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We report Molecular Dynamics (MD) simulations of a generic hydrophobic nanopore connect- 
ing two reservoirs which are initially at different Na + concentrations, as in a biological cell. The 
nanopore is impermeable to water under equilibrium conditions, but the strong electric field caused 
by the ionic concentration gradient drives water molecules in. The density and structure of water in 
the pore are highly field dependent. In a typical simulation run, we observe a succession of cation 
passages through the pore, characterized by approximately bulk mobility. These ion passages reduce 
the electric field, until the pore empties of water and closes to further ion transport, thus providing 
a possible mechanism for biological ion channel gating. 

PACS numbers: 61.20.Ja, 68.08.Bc, 87.16.Ac, 87.16.Uv 



Water and ion permeation of nanopores is a key issue 
for biological membrane-spanning ion channels and aqua- 
porins, as well as for materials like zeolites, gels and car- 
bon nanotubes. Recent simulations report intermittent 
filling of hydrophobic nanopores by water under equilib- 
rium conditions 0.0- Q However, imbalances in ion con- 
centrations in the inside and outside of cell membranes 
create strong electric fields |4|. Experiments on water 
near interfaces show that strong fields can induce con- 
siderable electrostriction of water The nonequilib- 
rium behavior of confined water and ions in strong fields 
should therefore be very important for ion permeation 
and ion channel function. In particular, biological ion 
channels "open" and "close" to ion transport in response 
to changes in the electric field across the membrane. This 
behavior, know as voltage gating, is crucial to their func- 
tion, but its mechanism is so far not well understood . 
Related approaches @ to ion transport include the appli- 
cation of a uniform external electric field Ljll£l : more spe- 
cific models of particular proteins 0, ITfl llll Il2| , Brow- 
nian Dynamics ,g| and continuum theories |l3| . Here we 
present the results of MD simulations in which a strong 
electric field across the pore is explicitly created by an 
ionic charge imbalance, as in a cell. We follow the relax- 
ation of this nonequilibrium system to equilibrium. 

Our generic model ion channel consists of a cylindri- 
cal hydrophobic pore of length L p = 16A and radius 
R p = 5 to 7A, through a membrane slab which separates 
two reservoirs containing water and Na + and Cl _ ions, as 
shown in Fig. 1. One reservoir has initial concentrations 
c Na + w 0.9M (12 cations) and c C \- « 0.6M (8 anions), 
while for the second reservoir, c Na + « 0.3M (4 cations) 
and cq\- ss 0.6M (8 anions). This imbalance of charge 
generates an average electric field of 0.37V/A across the 
membrane. These ion concentrations and electric field 
are typically five times larger than under "normal" phys- 
iological conditions, but they could be achieved in the 
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course of a rare, large fluctuation at the pore entrance. 
The enhanced ion concentrations in the initial state were 
chosen to improve the signal-to-noise ratio in the simula- 
tion and to allow the detection of novel transport mech- 
anisms. The chosen pore dimensions are comparable to 
those of the selectivity filter of a K + channel [3, 0] • 

The simulation cell contains two pores in sequence 
along the z-axis, one of which is shown in Fig.l. The 
reservoir to the right of this pore thus forms the left 
reservoir for the other pore. Due to periodic boundary 
conditions, the right reservoir of the latter is also the left 
reservoir of the first. In this arrangement, the simplest 
to allow the use of full three-dimensional periodic bound- 
aries, the flows of ions in response to the concentration 
gradient will be anti-parallel in the two channels. The 
relaxation towards equilibrium, where the two reservoirs 
are individually electroneutral, will thus involve an indi- 
rect coupling between the two pores. 

The water molecules are modelled by the SPC/E po- 
tential which consists of an O atom, carrying an 
electric charge q = — 0.8476e, and two H atoms with 
q = 0.4238e. The O-atoms on different water molecules 
interact via a Lennard- Jones (LJ) potential with pa- 
rameters e = 0.6502kJmol- 1 and a = 3.169A . The 
model is rig id, with OH bond length lA and HOH 
angle 109.5°. The Na + ions have LJ parameters e = 
0.3592kJmor 1 , a = 2.73A and q = +e and for the CI", 
e = 0.1686kJmor\ a = 4.86A and ?=-e 0. 

Ions and water O-atoms interact with the confining 
pore and membrane surfaces by a potential of the LJ 
form V = e'[(cr'/r) 12 — (er'/r) 6 ], where r is the orthogonal 
distance from the nearest confining surface. The poten- 
tial parameters are e' = 1.0211kJmor 1 and a' = 0.83A . 
If R p is the geometric radius of the cylindrical pore, one 
may conveniently define an effective radius R by the ra- 
dial distance from the cylinder axis at which the interac- 
tion energy of a water O-atom with the confining surface 
is fceT, leading at room temperature to R « R p — 2A ; 
similarly the effective length of the pore is L ss L p + 4A . 
Ion-water, ion-ion, water-surface and ion-surface cross 
terms are defined using the usual Lorentz-Berthelot com- 
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billing rules. Polarizability of the membrane [l7| an d of 
the water molecules and ions is neglected. 

The total simulation cell including both of the channels 
is of dimensions l x = l y = 23.5 ± 0.3A and l z — 112.9 ± 
1.5A and contains 1374 water molecules, 16 Na + and 16 
Cl~ ions. Molecular Dynamics simulations were carried 
out with the DLPOLY2 package using the Verlet 
algorithm |l9l . |2fJ with a timestep of 2fs. The pressure 
was maintained at P — lbar and the temperature at 
T = 300K using a Berendsen barostat and thermostat 
|2l|. Electrostatic interactions were calculated using the 
particle- mesh Ewald method |22^ . 






FIG. 2: Electric field in and around the pore, (a), 
Magnitude and direction of the electric field, depicted 
by vectors on a two dimensional rectangular grid in 
the x-z-plane. The contour of the membrane pore 
is sketched as black solid line. (b), Averaged z- 
component E z (z) of the electric field inside the channel. 



FIG. 1: Simulation snapshots. Molecular configura- 
tions (a) before a cation (dark grey spheres) perme- 
ates the channel and (b), 10 ps later, while it goes 
through. Anions are shown as light grey spheres. Only 
half of the periodically repeated simulation cell is shown. 



Water permeation of the pore is strongly affected by 
the electric field. The effective channel radius chosen 
for most of the simulations (i? = 3A) is such that under 
equilibrium conditions (i.e. with equal numbers of anions 
and cations on both sides and hence no electric field), 
the channel is empty of water and ions 4]. However, the 
ionic charge imbalance across the the membrane causes 
the pore to fill spontaneously with water. The electric 
field throughout the system was monitored by measur- 
ing the electrostatic force on phantom test particles on 
a three dimensional grid 0. Fig. 2a shows the average 
local electric field around one pore before the first ion 
moves through a channel. It is nearly zero in the reser- 
voirs. Inside the pore the field is very strong (~ 0.37V/A) 
and has a small inward radial component. The profile of 
the z-component of the field E z is shown in Fig. 2b and 
is constant inside the pore. 

During the course of the simulation, a number of Na + 
ions move through the pore. Each of these events changes 
the reservoir charge imbalance and reduces the electric 



field in the pore. This has a dramatic effect on the be- 
havior of the water, as shown in Fig. 3 in which the num- 
ber riH 2 o of water molecules inside one pore is plotted as 
a function of time for a typical simulation run. 

Initially, the water in the pore undergoes strong elec- 
trostriction, comparable to experimental observations 0, 
with an average density 

P ~ ?iH 2 o/(7ri? 2 £), (1) 

twice as large as that of bulk water in equilibrium. If 
we assume bulk density of water po inside a channel of 
radius R and length L we expect an average number ribuik 
of molecules inside the channel with 

^buik = po(nR 2 L). (2) 

"■bulk is indicated with an arrow on the right-hand side 
in Fig. 3. At each ion crossing, the number «h 2 o of 
water molecules inside the pore drops but is still larger 
than ribuik- In this particular simulation run shown in 
Fig. 3 three cations went through this pore. After the 
third ion crosses (t w 2.25ns), the electric field is no 
longer strong enough to sustain channel filling and the 
pore spontaneously empties of water, thereby becoming 
impermeable to further ion transport. However the other 
pore in the simulation cell remains filled and the final 
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FIG. 3: Occupancy and structure of water inside the 
pore. The number of water molecules inside the channel 
n n 2 is plotted as a function of time t. The shades of 
grey code the average magnitude of the electric field in- 
side the pore: from black to light grey: E ps 0.37y/A, 
E ps 0.25F/A, E Pa 0.15F/A , E ps 0.10V/A. The in- 
set shows the corresponding radial density profiles of the 
water molecules inside the channel averaged over periods 
of constant field, po is the bulk density of water in equi- 
librium, nbuik (eq. (2)) is the expected number of wa- 
ter molecules inside the pore if bulk density is assumed. 



ion crossing eventually occurs through this pore which 
after that also empties of water. Finally equilibrium is 
restored to the system and both channels are empty of 
water. On repeating the simulation 5 times we observe 
that the closing of one pore after the third ion passed 
through it occurs in all runs. 

The structure of water in the filled pore is strongly 
affected by the field, as shown in the inset of Fig. 3. 
Before the first ion crossing, water forms clear layers near 
the pore wall and along the z-axis. The central layer 
disappears after the first ion crossing and the outer layer 
becomes less well-defined. After a further ion crossing, 
water is more-or-less evenly distributed. 

Ion transport through the pore is found to occur es- 
sentially at constant velocity. Fig. 1 shows snapshots 
from a typical simulation run just before (a) and while 
(b) a sodium ion passes through the channel. Within the 
reservoirs the anions and cations diffuse among the wa- 
ter molecules. When a cation in the Na + -rich reservoir 
comes close to the channel entrance, it experiences the 
strong axial field shown in Fig. 2 and is dragged into 
the channel. Analysis of 15 simulation runs, with the 
same initial charge imbalance but different initial config- 
urations, shows that once the first ion enters the channel, 
it moves with a constant velocity which is approximately 
the same in all runs, and then reverts to diffusive motion 
at the other end of the channel. 

Fig. 4 shows typical cation positions along the z-axis, 
as a function of time, for the first, second and third ion 



FIG. 4: Cation trajectories inside the pore and in its 
vicinity. The z-coordinate of the cations is plotted versus 
time for three successive permeations in one typical simu- 
lation run. The channel is located between z = — 8A and 
z = 8A marked by the two long dashes lines, t = defines 
the time at which the ions are located at the entrance of 
the channel. The shades of grey code the average magni- 
tude of electric field experienced by the ions, as in Fig. 3 



crossings, as shown in Fig. 3. The second ion also tra- 
verses the channel at constant, somewhat reduced ve- 
locity, although it appears to pause for approximately 
lOps at the channel entrance, perhaps in order to shed 
its bulk-like solvation shell. We observe that this "paus- 
ing time" is rather widely distributed between simulation 
runs. The cation mobility defined by v = /i+eE, can 
be calculated from the slopes of the trajectories in Fig. 
4, together with the measured electric fields, as in Fig. 2. 
The resulting values are /i+ ps 4.5 x 10 11 skg -1 for the 
first ion, pi + ps 3.8 x 10 11 skg"" 1 for the second ion, and 
fi+ ps 2.4 x 10 11 s kg -1 for the third ion. These values are 
close to the value of /i+ ps 2.3 x 10 11 s kg -1 obtained from 
the self-diffusion constant in the reservoir D + , using Ein- 
stein's relation /i + = D^/ksT, but seem to increase with 
the magnitude of the electric field inside the pore. This 
enhancement of the mobility correlates with the change 
of structure of the water inside the channel, shown in 
the inset of Fig. 3. The tetrahedral hydrogen bond net- 
work which water forms under equilibrium conditions is 
disrupted inside the pore under high electric fields. 

Simulations of wider (R — 5A) pores and with different 
lengths give qualitatively the same results. The critical 
electric field for water permeation is, however, sensitive 
to the pore radius and length. This suggests that voltage- 
dependent gating in ion channels, if it were to occur by 
changes in water permeation of a hydrophobic section of 
the pore 0, H , might be strongly dependent on channel 
geometry. 

The key finding which emerges from our simulations 
is the strong correlation between water and ion behavior 
under non-equilibrium conditions. Ionic charge imbal- 
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ance across the membrane induces water permeation of 
the hydrophobic pore and thus makes it permeable to 
ions. This suggests that voltage gating of ion channels 
may be linked to the coupling between water and ion per- 
meation in pores far from equilibrium. The structure and 
density of water in the pore is dramatically affected by 
the strong electric field. The passage of a cation through 
the channel causes an abrupt jump of the electric field, 
and an ensuing jump in the number of water molecules 
inside the pore. Ion passage through the pore occurs at 
constant velocity and with a mobility coefficient similar 



to that of the bulk solution at equilibrium. 
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